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ABSTRACT 


A  non-linear  thick  composite  shell  theory  is  presented 
in  which  the  t h rough- t he - t h i c kne s s  displacements  are  model¬ 
led  using  a  variation  of  a  cubic  spline.  The  theory  is 
developed  by  considering  the  Lagrangian  strains  in 
conjunction  with  the  2nd  Piola-Kirchhof f  stress.  This 
formulation  leads  to  a  theory  which  encompasses  large 
displacements  with  moderately  large  rotations  but  is 
restricted  to  small  strains.  The  imposition  of  the  cubic 
distribution  t hr ough - t he - t h i c kne s s  insures  that  the 
compatabi 1 i ty  of  the  displacements  and  their  first  and 
second  derivatives  and  thus  the  shear  strains  are  maintained 
from  lamina  to  lamina.  The  cubic  distribution  is  seen  as  a 
higher  order  approximation  than  has  been  previously 
employed,  but  because  of  the  nature  of  the  spline,  the 
theory  is  less  cumbersome  and  more  easily  implemented  than 
the  parabolic  theory.  In  addition,  there  is  no  introduction 
of  additional  degrees  of  freedom  with  the  cubic  theory.  A 
family  of  2-D  isoparametric  elements  is  employed  in 
conjunction  with  the  theory  to  solve  a  class  of  3-D  thick 
plate  problems.  Results  are  presented  showing  comparisons 
which  are  in  good  agreement  with  previous  work. 


CHAPTER  I 


INTRODUCTION 


BACKGROUND : 

In  the  last  fifteen  to  twenty  years,  composites,  espe¬ 
cially  fiber  reinforced  laminates  have  found  increasing 
application  in  aerospace  structures,  pressure  and  underwater 
vessels,  and  nuclear  reactor  structures.  The  advent  of 
advanced  f i be r - r e i n f o r c e d  composite  materials  has  been 
called  "the  biggest  technical  revolution  since  the  jet  en¬ 
gine"  (1)  and  has  led  to  increasing  use  of  laminated  plates 
and  she  1 1 s . 

Fiber-reinforced  composites  such  as  boron-epoxy  and 
gr aph i t e -e poxy  possess  a  number  of  desirable  features.  Two 
I  are  their  high  stiffness  to  weight  ratios  and  their  aniso¬ 

tropic  nature.  These  traits  when  exploited  by  the  designer 
provide  the  means  for  constructing  weight  efficient,  strong 
structures  which  respond  to  loads  in  a  manner  unlike  iso¬ 
tropic  materials.  A  good  example  of  such  an  application  is 
the  forward  swept  wing  fighter  currently  in  the  development 
stage.  The  inherent  wing  divergence  problem  of  such  a 
design  is  controlled  by  the  judicious  lay-up  of  composite 


1  am  i  nae 


In  order  that  full  advantage  be  taken  of  these  desir¬ 


able  traits  ,  it  is  necessary  to  optimize  the  design.  This 
optimization  process  requires  the  use  of  valid  theories  to 
predict  structural  response  to  a  variety  of  load  conditions. 

Many  of  the  early  theories  pertaining  to  the  analysis 
of  laminated  plates  and  shells  were  based  on  the  classical 
Kirchhoff  hypothesis  and  are  well  established  in  the  lit¬ 
erature  (2-6).  Exact  elasticity  solutions  for  some  part¬ 
icular  plate  bending  problems  have  been  obtained  by  Pagano 
(7-9)  and  Srinivas  and  Rao  (10).  Comparisons  of  results 
obtained  from  the  classical  theory  with  the  exact  elasticity 
solutions  indicate  the  necessity  of  considering  transverse 
shear  deformations  in  the  analyses  of  laminated 
s t r uc t u r e s . ( 1 1 )  This  necessity  is  implied  because  the 
distortion  of  the  deformed  normal  due  to  transverse  shear 
strain  is  dependent,  not  only  on  the  laminate  thickness,  but 
also  on  the  orientation  and  degree  of  orthotropy  of  the  ind¬ 
ividual  layers  (11). 

One  of  the  earliest  attempts  to  include  the  shear 
strain  effects  for  anisotropic  laminated  plates  was  made  by 
Ambartsumyan  (12).  He  assumed  that  the  thickness  shear 
stress  was  distributed  as  a  single  continuous  parabolic 
function  through  the  entire  laminate  ,  Just  as  in  a  homo¬ 
geneous  shell.  It  was  shown  by  Hsu  and  Wang  (13)  that 


Ambartsumyan's  assumption  was  not  valid  for  a  laminated 
shell  (since  it  does  not  permit  satisfaction  of  interlaminar 
compatibility)  unless  the  in-plane  Poisson's  ratios  are 
Identical  in  all  of  the  layers.  Yang,  Norris  and  Stavsky 
also  developed  a  theory  for  heterogeneous  plates  (14)  in 
which  they  assumed  that  the  in-plane  displacements  varied 
linearly  with  the  thickness  coordinate.  They  had  to  use 
correction  factors  to  satisfy  the  boundary  conditions. 
Whitney  applied  the  work  of  Ambartsumyan  and  Yang  to 
specific  laminated  plate  problems  (15-17).  He  also  present¬ 
ed  a  theory  for  the  analysis  of  cylindrical  shells  which 
included  both  transverse  shear  and  transverse  normal  strain 
(18).  The  model  used  was  a  direct  extension  of  the 
Hildebrand  et.al.  theory  (19,20)  for  homogeneous  isotropic 
shells.  Further  extension  of  the  Hildebrand  theory  was  per¬ 
formed  by  Lo  et.al.  (21,22).  In  the  Lo  theory,  in-plane  and 
transverse  displacements  were  assumed  to  be  cubic  and 
quadratic  functions,  respectively,  of  the  through-the- 
thickness  coordinate. 

The  next  logical  step  following  the  work  of  Yang, 
Norris,  and  Stavsky,  was  to  assume  that  the  displacements 
varied  linearly  within  each  layer  of  the  laminate.  Such 
work  was  performed  by  Grot  (23)  and  by  Srinivas  (24).  An 
approach  similar  to  that  developed  by  Srinivas  was  used  by 
Epstein  to  formulate  the  equations  for  the  nonlinear  anal- 


y  6 1 s  of  multilayered  shells  (25). 

The  finite  element  method  of  analysis  for  laminated 
plate  bending  problems  was  presented  by  Pryor  et.al.  (26). 
This  formulation  assumed  a  uniform  shear  strain  angle 
through  the  thickness  of  the  plate  and  neglected  local 
effects  such  as  the  state  of  stress  and  deformation  at  the 
layer  interfaces.  Perhaps  the  first  attempt  at  considering 
these  local  effects  was  undertaken  by  Mawenya  et.al.  (26), 
wherein  Ahmad's  s u pe r pa r ame t r i c  quadratic  shell  element  (27) 
was  used  with  independent  normal  rotations  within  each 
layer . 

Panda  et.al.  (28)  combined  the  ideas  of  Pryor  and 
Mawenya  by  again  using  the  Ahmad  element  as  did  Mawenya  but 
maintaining  a  uniform  shear  strain  angle  as  did  Pryor.  They 
later  applied  their  theory  to  shell  structures  (29),  using 
a  doubly  curved  Ahmad  element. 

With  the  exception  of  Epstein's  (25)  work,  each  of  the 
preceding  developments  was  for  small  displacement  linear 
theory.  The  work  done  by  Witt  (30)  appears  to  have  been  the 
first  attempt  at  satisfying  the  continuity  of  strains  at 
layer  boundaries  while  using  a  nonlinear  large  displacement 
theory.  He  assumed  that  the  normal  rotations  varied  linearly 
within  each  layer.  Thus,  the  displacements  varied 
quadrat i ca 1 1 y  .  This  was  a  direct  extension  of  the  work  by 


Mawenya  (26),  but  employed  a  simpler  triangular  element  and 
a  fully  nonlinear  formulation. 

In  more  recent  years  a  number  of  investigators  have 
studied  various  methods  of  including  the  transverse  shear 
strain  effects.  Pack  and  Man del  (31)  employed  a  concept 
whereby  they  derived  a  "connector"  element  which  essentially 
linked  2-D  elements  in  an  attempt  to  effectively  model  a  3-D 
structure.  Owen  and  Figueiras  (32)  employed  a  degenerate  9- 
noded  3-D  continuum  element  (27)  in  which  they  assumed  a 
constant  transverse  shear  strain.  A  correction  shear  factor 
derived  in  cylindrical  bending  was  used  to  approximate  the 
real  shear  strain  energy  component. 

The  extension  of  Panda's  work  to  the  nonlinear  theory 
was  performed  by  Chang  et.al.  (33).  He  used  the  doubly 
curved  Ahmad  element  in  conjunction  with  the  updated 
Lagrangian  description  of  the  virtual  work  equation. 

The  use  of  the  so-called  hybrid  or  mixed  finite  element 
method  has  been  explored  by  a  number  of  researchers.  These 
applications  range  from  linear  plate  theory  (34-38)  to  non¬ 
linear  shell  theory  employed  by  Noor  et.al. (39).  Noor  has 
extended  his  work  to  collapse  anlysis  of  shells  (40). 

During  the  course  of  this  work,  a  number  of  other 


authors  have  published  papers  relevant  to  the  area  of 
research.  Reddy  and  Creamer's  (41)  work,  was  an  analytical 
approximation  approach  in  which  they  expanded  the  In-plane 
displacements  as  cubic  functions  of  the  thickness  coordinate 
and  assumed  the  transverse  deflections  to  be  quadratic 
through  the  thickness.  Solutions  to  the  simply  supported 
flat  plate  equations  were  obtained  using  Navier's  method. 

Epstein  and  Glockner  (42)  specialized  their  previous 
work  (25)  to  a  linear  analysis  of  deep  and  multilayered 
beams  —  solving  the  problem  using  a  finite  difference  method. 
Epstein  and  Huttelmaier  (43)  further  specialized  the  Epstein 
theory  to  thick  multilayered  plates.  In  this  work  they 
employed  a  4-noded  isoparametric  finite  element  to  solve  a 
linear  plate  bending  problem. 

The  use  of  3-D  elements  to  model  each  individual  lamina 
has  been  done  for  a  number  of  years.  The  most  recent  works 
have  been  those  of  Lee  (44)  and  Kuppusamy  and  Reddy  (45). 
In  these  works,  the  8-noded  trilinear  elements  were  used  for 
linear  as  well  as  nonlinear  analyses.  The  use  of  these 
elements  is  restricted  to  a  class  of  small  problems  due  to 
the  large  number  of  degrees  of  freedom  needed  to  effectively 
model  a  multilayered  plate. 


PROBLEM  STATEMENT: 


In  the  present  investigation,  a  family  of  finite  ele¬ 
ments  is  constructed  which  accounts  for  the  following  fea¬ 
tures:  (most  of  which  have  not  been  effectively  combined  in 
previous  works) 

1.  Nonlinear  geometrical  effects  related  to 
large  displacements  and  moderately  large 

r  o  t  a  t i on  s . 

2.  Variations  of  transverse  shear  strain  through 
the  thickness. 

3.  Heterogeneity  through  the  thickness. 

The  transverse  shear  strain  effects  are  included  by 
assuming  that  there  is  a  linear  distribution  of  transverse 
curvature  of  each  lamina,  thus  implying  a  parabolic  distri¬ 
bution  for  the  slope  and  a  cubic  distribution  for  displace¬ 
ments  through  the  thickness.  This  choice  of  function  is  a 
natural  extension  of  the  previous  work  (in  particular  that 
of  Witt). 

T'te  study  is  limited  to  a  flat  plate  having  individual 
laminae  of  equal  thickness  but  varying  material  properties 
and  fiber  orientations. 


The  derivations  that  lead  to  the  finite  element  stiff¬ 


ness  formulation  are  discussed  in  Chapter  II.  A  computer 
program  based  upon  the  results  for  a  flat  plate  is  discussed 
in  Chapter  III.  The  results  for  several  problems  investiga¬ 
ting  the  applicability  and  range  of  the  elements  are 
discussed  in  Chapter  IV.  A  number  of  conclusions  are 
presented  in  Chapter  IV  as  well.  These  conclusions  are 
reviewed  and  additional  conclusions  are  presented  in 
Chapter  V. 


a 


CHAPTER  II 


THEORY 


GENERAL  SHELL  EQUATIONS: 

In  this  section,  the  objective  is  to  present  the  dev¬ 
elopment  of  the  general  shell  equations  from  a  kinematic 
viewpoint,  lay  out  the  technique  used  to  maintain  strain 
continuity  at  the  interlaminar  boundaries,  and  cast  the 
theory  in  finite  element  terms. 

The  derivations  of  the  general  shell  equations  from  the 
Lagrangian  standpoint  is  a  straight-forward  process  which 
has  been  carried  out  in  a  number  of  previous  works  (for 
example  Saada  (46))  but  is  presented  here  for  completeness. 
Also,  the  assumptions  which  are  made  in  the  derivation 
differ  from  those  of  Saada  in  that  a  line  element  which  is 
normal  to  the  shell  reference  surface  is  not  allowed  to  vary 
in  length.  The  equations  which  are  finally  presented  will 
permit  analysis  of  a  class  of  problems  which  are 
characterized  by  large  displacements,  moderately  large 
rotations,  and  small  strains.  One  of  the  main  reasons  for 
choosing  the  Lagrangian  strain  for  this  work  was  that  the 
material  properties  could  be  referenced  to  the  undeformed 
con f i gu r a t i on  and  not  varied  in  the  analysis.  Material 


nonlinearity  is  not  considered  in  this  work,  but  could  be 
easily  employed  in  future  work. 


Referring  to  Fig.  1,  consider  the  displacement  f 
de  fined  by  : 


*  =  i  '  3 


and  the  base  vectors  in  the  undeformed  space  as: 


*  -  7% 

The  base  vectors  for  the  deformed  space  are  given  by 

?  -  4  - « '  4 ) 


The  metrics  for  the  two  spaces  become: 


Aj  =  A-  •  4 


for  the  undeformed  and 


r*  st  * 


for  the  deformed  space. 


The  definition  of  the  Lagrangian  strain 
(which  is  a  strain  measure  with  reference  to  the 
configuration)  is  given  by  Saada  (46)  as: 


tensor,  j  ■■  , 

vJ 


undef  ormed 


( 


With  proper  substitution  of  the  metric  tensor  and  simple 
algebraic  manipulation,  one  obtains  the  expression 


The  unit  vectors  written  in  terms  of  base  vectors  are.' 

(no  surnnat i on ) 

Z  _  Z_ 

^  /  At  I  /  All 

thus  (adopting  a  summation  convention) 

A  =  hi  e, 

The  displacement  vector  written  in  terms  of  these  unit 


vectors  becomes 


Substitution  of  (9)  and  (10)  into  the  definition  of  the 


Lagrangian  strain  (7)  yields 


2  {'«.•'  =  kei  ■  ‘fyy.UJk)  J-hj-ej-  %j.  (u*e.J)  + 

U  J 


From  the  study  of  the  theory  of  surfaces,  one  sees  that 
any  point  on  a  shell  can  be  located  by  means  of  three 
parameters.  two  of  which  vary  along  a  reference  surface 
while  the  third  varies  along  the  reference  surface  normal. 
The  lines  of  curvature  of  the  reference  surface  are  chosen 
as  parametric  curves.  These  lines  together  with  the  normal 
form  an  orthogonal  system  of  reference.  An  arbitrary  point 
in  the  shell  is  located  by  means  of  the  position  vector 

where  IT  is  the  position  vector  of  a  corresponding  point  on 
the  reference  surface  and  y3  is  the  distance  of  the 
arbitrary  point  from  the  reference  surface  measured  along 
the  unit  normal,  t (See  Fig  la.)  The  magnitude  of  an 
element  of  length  is  given  by: 

dA.  =  (J$  ■ 

Following  manipulation*^  is  written  as: 


1  4 


(14) 


25  -4,2,4',  + $>  Jys 

where4A  *  =  X/c<  (the  base  vectors  for  the  ref.  surface) 

From  Saada  (46),  one  sees  that  the  derivatives  of  the 


unit  vectors  are: 


Je,  - 

4f' 


—  s 

ft  J 


(15) 


c)  £  i  _ 

Sid 

*1 

(16) 

it 

W  1 

^  -T’ 

^'i2-  c, 

(17) 

il 

—  z-  Jo"  _ 

JU.  3 

f  ^ 
c —  1 

(18) 

U  Cw 

~  \a 

It 

^  1  /O 

-e,  ' 

(19) 

^<r3  _ 

«*  5 

A. 

(20) 

re  Rj  and 

represent  the 

principal  radii  of 

curvature 

and  derivatives  with  respect  to  are  assumed  to  be  zero. 


1  5 


Substitution  of  the  appropriate  derivatives  into  eqn. 


(14)  gives: 

Ay  -  (a,  Ayt  +j/3 


so 


t  hat 


Since  by  definition  (Saada  (46)): 


it  is  seen  that 

h,=  a,(i  oy*,) 

K  -  i  *-  y/L) 

h  3  -  / 


Substitution  of  (24-26)  and  (15-20)  into  the  definition  for 


Lagranglan  strain  yields  expressions  for 


%  ' 


The  full 


derivation  is  given  in  Appendix  I. 


(21  ) 


(  22  ) 


(23) 


(24) 


(25) 


(26) 
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The  displacement  vector  is  divided  into  two  parts,  (see 

Fig.  la)  one  part  Is  a  ;  a  function  only  of  the  surface 

coordinates  and  associated  with  the  deformation  of  the 

A 

reference  surface.  The  other  part  i s  ;  a  function  of 

all  three  coordinates  and  associated  with  the  deformation 
through  the  thickness.  The  components  of  u  are  written 
a  s  :  (  3  0  ) 


Hi 


u.t 


<y,>Ki  y») 


««.  =  ft 

Key., >o 

so  that 

3ii  r  4>0  +liVe,)[(u,  57  +  Ui  jrji  * 

’k^LcX'+Z'),,  +-  + 

[  -  (u>  +u)  + 


(27) 


(28) 


(29) 


(30) 


4- 


A  comparison  with  Witt's  work  (30)  reveals  that 
although  the  method  of  derivation  is  different  here,  the 
same  results  are  obtained. 


Conversion  from  tensorial  to  physical  strain  is  performed 
so  that  material  properties  can  be  readily  found  and  used. 
This  is  accomplished  from:  (46) 


So  that 


e 


i 


IM 


&  -  M 1  *} 


(Note  that  u,v,w  are  adopted  here  to  represent  the  displace 
ments  u^,  u^,  and  u^.  This  is  done  to  avoid  confusion  i 

subsequent  derivations.) 

The  [L]  operator  matrix  is  further  reduced  to  it 
component  parts: 

[l]-[lo]+[li]  +[l2] 

where  LO  contains  only  linear  expressions,  LI  contains 
nonlinear  expressions  involving  only  the  reference  surface 
displacements,  L2  contains  nonlinear  expressions  involving 
only  the  through  the  thickness  displacements. 

In  order  that  the  L's  could  be  separated  out,  it  was 
necessary  to  expand  out  the  expressions  for  the  strains 
The  full  expansion  was  carried  out  and  the  resulting  (L 


components  are  found  in  Appendix  II.  The  expressions  found 
in  Appendix  II  are  exactly  the  same  as  Witt's  (30)  but 
differ  only  in  form. 


CUBIC  SPLINE  FORMULATION: 

In  Witt's  work  (30),  he  considered  the  use  of  a 
parabolic  assumption  for  the  distribution  of  u  through  the 
thickness  of  the  shell.  This  was  done  by  assuming  that  the 
rotations  were  distributed  linearly  through  each  lamina. 
The  parabolic  nature,  of  u  was  then  obtained  by  integrating 
the  rotations  and  satisfying  continuity  of  u  and  the 
rotations  at  the  interlaminar  boundaries.  With  this 
assumption,  he  was  able  to  show  that  strain  compatibility 
could  be  maintained  at  the  lamina  interfaces.  Accurate 
results  were  obtained  when  each  lamina  was  modeled  as  two 
or  more  "pseudolayers".  This  modeling  resulted  in  an 
increase  in  the  number  of  degrees  of  freedom  and  thus 
increased  the  computer  time  and  storage  needed  for 
so  1 u  t i on . 

In  this  work,  an  attempt  was  made  to  increase  the  order 
of  the  assumption  for  the  variation  of  u  through  the 
thickness  from  parabolic  to  cubic.  This  was  done  to  more 
accurately  represent  the  real  displacements  without 
resorting  to  the  inclusion  of  psoudolayers  with  their 
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inherent 


increase  in  degrees  of  freedom. 


A  natural 


extension  from  Witt’s  work  was  to  begin  with  a  linear 
assumption  for  the  second  partial  of  u  with  respect  to  the 
thickness  coordinate,  integrate  twice,  and  evaluate  the 
constants  of  integration  so  as  to  insure  strain  compatibil¬ 
ity  through  the  shell.  This  thinking  led  to  a  variation  of 
cubic  spline  theory  in  which  the  u's  were  represented  as 
functions  of  the  rotations.  The  following  section  shows  how 
this  thinking  was  implemented. 

.  A 

Concentrating  our  attention  on  the  u  component,  we 
consider  now  the  shell  of  concern  to  be  made  up  of  M  layers 
of  orthotropic  laminae  each  having  its  own  set  of  material 
properties  which  may  or  may  not  be  unique  within  the  shell 
(see  Fig*  2).  We  are  concerned  that  as  we  observe  the 
interlaminar  boundaries  that  the  strains  be  continuous  from 
one  lamina  to  another. 


Le  t 

the 

bottom 

of  the 

jth  layer  be  given 

the  coordinate 

z^  j  and 

the 

top 

V 

whe  r  e 

z  “  y  •}  •  At 

•7Z 

; 

/« 

-  Kj., 

( 

and  a  t 

• 

a  u 

-  y 

c 
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If  the  thickness  of  the  jth  layer  is  denoted  by 


Zj  ~  -l 


and  within  that  layer,  a  linear  assumption  is  made  for  the 
distribution  of  the  second  partial  of  u  with  respect  to  z, 
then 


Jji  =  1/  ,  n 

CHT  'V'  tfy 


( 


This  assumption  is  recognized  as  the  well-known 
starting  point  for  the  development  of  cubic  spline  theory, 


wh ich  is  usually 

used 

in  connecting  graphi 

cal 

points  wi 

t  h  a 

smoot  h  curve. 

A  number  of 

references 

are 

replete 

with 

documentation  of 

this 

theory . 

One  of  the 

best 

t  r eatment  s  is 

found  in  reference  (47). 


If  we  integrate  eqn  (50)  and  evaluate  the  constants  of 

integration,  we  obtain  the  equations:  (see  reference  47) 


-l 


<w> 


•  *$■ 


ztj 


A  A 

«,•  -K,w 


-  i 


t)' 


V 


(51  ) 


(52) 


A  A 

where  is  the  value  of  u(z)  at  z  =  z^. 


The  integration  process  yields  the  desired  parabolic 
distribution  for  the  first  derivative  and  a  cubic  distribu¬ 
tion  for  the  displacement  function. 


At  this  juncture,  we  impose  the  continuity  requirement 


^  l*j 


2 


(which  means  that  the  slopes  are  the  same  as  one  approaches 
the  interlaminar  boundary  from  top  or  bottom). 


and  obtain  the  condition: 


A)S<  ^ 


-f%  -«j.i)/b-3 


ti  * 


where : 


j  .  s  JilL  • 
J  ij  *  tl+i  ) 


\  A  ■-  H  (s‘ 


Equations  (52)  and  (54)  could  be  used  to  develope  an 

A 

interpolation  relation  for  u.  Such  an  interpolation  would 

A  A1* 

involve  the  second  partials  of  u  with  respect  to  z  \-t--t)  a  t 

d* 

each  interface  and  would  represent  a  sound  theory  for 
implementation  in  finite  elements.  In  other  words  looking 
ahead,  some  of  the  nodal  degrees  of  freedom  for  the  element 
would  be  the  curvatures  at  the  interfaces.  Since  Witt  used 
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the  rotations, 


and  since  rotations  are  easily  interpreted 


from  a  physical  standpoint,  the  foregoing  development  is 


recast  in 

terms 

of  the 

rotations  or  slopes  as  follows 

2  •= 

ii-' 

) 

—  -  / 

ai 

^  - 

sj- 

. 

iA  _  J 

) 

Equations  (51,52,54)  become: 


Equations  (54  and 
systems  of  M- 1  equations 


59)  represent 
in  the  M+ 1  unknowns 


underdetermined 

A 


One  of  the  unknowns  is  eliminated  by  choosing  the 


(56) 


(57) 


(58) 


(59) 


reference  surface  to  be  at  the  inner  surface  of  the  shell 


and  thus  u^=0. 

The  remaining  necessary  condition  to  make  the  system 
tractable  must  be  imposed  at  the  outer  surface  of  the 


shell.  This  condition  was  dealt  with  by  assuming  an  imagi¬ 
nary  layer  to  exist  at  the  upper  surface.  By  allowing  no 
changes  in  rotation  and  no  deformations  to  occur  in  the 
imaginary  layer,  the  u n d e r d e t e r m  i  n e d  system  becomes  tract¬ 
able.  A  simple  2-laver  example  of  this  procedure  is 
presented  to  make  the  process  easily  understood.  This  2- 
layer  case  is  then  generalized  to  an  M-layer  case. 

As  an  illustration  of  how  the  final  condition  is 
imposed  ,  a  two-ply  laminated  plate  as  shown  in  Fig.  (  3  ) 
is  considered.  Note  that  if  a  third  (non-existent)  layer  is 
imagined,  the  following  set  of  equations  is  obtained  by 
direct  substitution  into  eqn  (59): 


o  r 

+  ^  (u,  -U.) 

/D  :  %  («*,-%■) 


(60) 


o\o\o\o 


o 


written  in  vector  form 


A 

Next,  the  reference  surface  boundary  condition  Ug*=0  is 
imposed.  The  boundary  condition  at  the  outer  surface  is 
expressed  as  and  U-u-f  It  .  These  conditions 
account  for  the  non-existent  layer  and  show  that  a  line 
extending  into  the  imaginary  layer  is  not  deformed,  but 
follows  the  natural  shape  imposed  by  the  existent  layers. 
In  spline  theory,  this  is  analogous  to  the  natural  cubic 
spline  condition.  ( 4 7  ) 


The  set  of  equations  thus  becomes 


141  7.  ,  %  ° 


L.o  i  ^  L-i  n  u,}  <■ 

Which  is  tractable  since  there  are  2  equations  and  2 


unknown  s  . 


Switching  to  a  more  general  matrix  notation  yields  the 
system: 


L*Wl  ‘  %  M/s] 

Solution  of  the  system  for  <j[>  gives: 

41=  %  [4'ls]  41  --  %  Lci  Ml 

where 

m'  [ ; ;']  ,  w -  m'k 

Written  out,  the  solution  is: 


(64; 


(65] 


(66) 


(67 


These  values  for  <u>  are  substituted  into  eqn  (58)  to 
yield  the  cubic  spline  approximation  for  the  displacements 
in  the  layers  in  terms  of  the  interlaminar  rotations  <^>. 


If  we  define  s(z)  as  the  cubic  approximation  for 


A 

u  (  z  ) 


in  the  Jth  layer,  we  see  that  for  this  two  ply  case  the 


V/2--  i-  fair*')  **]  s  ? 

/  *ffc  ^  *  3^2  ^3 


///3  =■ 


fgr-j/)  [2(2,-*)  1  £ 

3** 


s&t^)z  ~  /4/  ^  +-  //*2  *fL 


z  where  a 

,  ^JA  +  te-z7)b(u-z)rt]„  \ 

™  I  3t"  "  3^=-  Z'J 


/  -  f  Ci-ih(ix-i)L r  +  (i 


(The  ' s  are  the  elements  of  the  (Cj  matrix  of  eqn  (67)) 

If  one  were  to  go  through  this  exercise  for  a  succes¬ 
sively  Increasing  number  of  plies,  he  would  see  a  recursion 
relation  emerge  and  would  find  that  in  general: 

a*  -^r-  ■/* 

C:  ,  '  +  £.  .  ( 

J_,‘  M*  J<" - W*  j 

whe  r  e  Xv  is  Kronecker's  delta  and  0^=0. 
vJ 

Using  the  expression  for  H  ,  one  may  write  the  inter¬ 
polating  spline  in  the  jth  layer  as: 

Af  +  f 


or  (adopting  a  summation  convention)  as: 


-  tijt  *  *  *  M  +  l) 


In  this  form,  the  H  can  be  thought  of  as  a  set  of  shape 

functions  in  the  z  direction.  It  should  be  noted  that  this 
derivation  for  H  is  original  work  and  represents  a 

variation  in  cubic  spline  theory  which  to  the  best  of  the 
author's  knowledge  has  not  been  previously  explored. 


v  component  of 

and  arguments  as 

for  H  ,  .  wh i c  h  are 

i  i 

FINITE  ELEMENT  APPLICATION: 


If  we  turn  our  attention  to  the 
displacement  and  make  the  same  assumptions 
we  did  for  u,  we  would  obtain  expressions 
identical  to  those  already  obtained. 


At  this  juncture  attention  is  concentrated  on  the  in¬ 
plane  directions.  Using  the  finite  element  shape  functions 
N^  for  the  in  plane  distributions,  we  write  (adopting  the 
summation  convention): 

-  A/*  $ L  *  ( 


whe  r  e 
node 


at  the 


the  partial  of  u  with  respect  to 
ith  interface.  Substituting  (75) 


z  for  the  kth 
into  (74): 


AU)j  n  tiji  hc-\ 


i  =-  1 ) 2 ,  •  •  •  M  +  | 

K.  =-  lf  2 1  .  .  .  Y\  neJes  ( 


Remembering  that 
in  the  jth  layer. 


j  - 


1  (  z  )  ^  is  the  cubic  approximation  for  u(z) 


we  recall  that  for  the  jth  layer: 


iu*)j  -  u  uc&j 

A 

Upon  substitution  for  u(z)  we  have: 


C 


11/ Sr),  ^  U  .• 


( 


The  forgoing  arguments  are  carried  out  on  the  v  dis¬ 
placements,  this  time  letting  be  the  partial  of  v 

with  respect  to  to  z  at  the  kth  node  and  ith  interface  so  that 

A 

one  obtains  a  similar  result  for  the  v  components  as  those 
for  the  u  components. 


C  \ 


Casting  all  of  the  preceding  development  in  finite 
element  terms  is  straight  forward.  We  choose  an  n-noded 
isoparametric  element  and  define  our  nodal  degrees  of  free¬ 
dom  for  the  kth  node  as: 


(79) 


This  choice  for  the  ordering  of  <a>  allows  for  variation  in 
the  element  type  that  is  very  efficient  in  terms  of  computer 
implementation . 


The  shape  functions  for  the  kth  node  are  then: 


tfj,  AV 
0  0 

o  0 

I  ^ 


pjn  A/«  0  O  •  •  •  0  O' 

o  a u  fad* 0 

0  O  O  •  *  •  O  tin. 


(80) 


» 


35 


where  the  H  's  are  given  by  eqn  (72)  and  the  N  's  are  the 
J  *  * 


in^plane  shape  functions. 


For  the  purposes  of  illustration  and  for  subsequent 
use,  the  isoparametric  quadrilateral  flat  plate  element  is 
considered.  The  shape  functions  are  given  by  (see  Fig.  4) 

A O  *$.)  ('  fp/^-  < 

where 

The  displacement  vector  is  thus  seen  to  be: 


At  this  point  in  the  development,  all  of  the  pieces  are 
available  for  derivation  of  the  finite  element  equations. 
This  derivation  is  readily  carried  out  by  application  of 
Hamilton's  Principle  which  states: 

l  J> 

/(  ( r-v)Ji  +  f  J"  w/  J*  =  o  < 


where  T  is  the  kinetic  energy,  V  is  the  strain  energy 
is  the  external  work. 


Looking  at  these  terms  in  detail,  it  is  seen  that 


where  <b>  is  a  vector  of  body  forces  and  <f)  is  a  vec 
surface  tractions. 


where  <S>  is  the  energy  complement  of  the  Lagrangian 
known  as  the  second  Piola-Kirchhof f  stress  vector. 
The  s t r e s s - s t r a  i  n  relation  is  given  by  <S>  =  [D]<e>. 


where  the  dot (  1  )  denotes  time  differentiation,  and 
matrix  of  mass  density  for  the  structure. 

Substitution  of  these  expressions  into  Hamilton 


ciple  yields: 


sCi 


4- 


where  t  i s  t ime  . 


[ 


s 


,  and  W 


(85) 

tor  of 

(86) 

strain 

(56) 

(87  ) 
p]  is  a 

Pr  in- 


(88) 


Substitution  for  <u>,<e> 


and  <S>  in  terms  of  (L],  [N], 


and  (a)  yields: 


a 


(89) 


Upon  taking  the  indicated  variations  and  reordering,  it  is 
seen  that: 

J  Mid  Col  MWm  I  - 

MTfC]  J  l'  <  ’°  > 

o  r 


MM  *  WW  =  [8]  +  ff] 


Equa  t i on 
structure  of 
limiting  our 
body  forces. 


(91)  represents  an  equation 
interest.  This  equation 
attention  to  static  problems 
Thu  s  : 


of  motion  for  a 
is  simplified  by 
while  neglecting 


/M  H 


(92) 
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The  definition  of  <a>  has  been  well  established 


1  n 


equation  (79).  The  formation  of  [K]  and  (F) in  eqn  (92)  are 
presented  to  provide  clarity  and  completeness. 


The  stiffness  matrix  ( K ]  is  given  by: 


£/v]  [l\  Co]  d/o] 


(93) 


Account  for  the  varying  material  properties  through  the 
thickness  as  well  as  the  changes  in  [N]  is  taken  by  forming 


the  sum  over  the  number  of  layers  as: 

[k]  =  £  f  ca/l  [l]tld]  o.]  [tJi  j  vi  / 

L  Jr,  J  v*/  J  J  J 

Since  the  operator  matrix  [L]  is  given 


(94) 

by 


[ L  ]  =  [ LO  )  +  [ L 1  ]  +  [ L2  ]  ,  further  factoring  of  [K]  is  accomplished 
as  (  K  ]  *  [  KO  ]  +  [  K  1  ]  +  [  K2  ]  where: 

[KOI  =  Z  [  [4  ^ CdJj  [loI  [»%■  J  V.  I  ( 


t  f  aaiuJj 

[>jJj  [uj  [DjjMWjj 


4 - 


[vjj  [U]T  [d]j  & J  La/Jj  J  %  I 


(95b) 


(95c 


This  way  of  writing  the  stiffness  matrix  lends  itself 
well  to  gaining  an  understanding  of  the  meaning  of  the 
various  terms.  The  [ K  0  ]  term  is  seen  to  contain  only 
linear  operators  [LO]  which  are  displacement  independent. 
Inclusion  of  only  this  term  in  the  computation  of  ( K  ]  is 
accounting  for  the  commonly  used  small  displacement,  small 
strain  components  in  the  matrix.  Further  inclusion  of  the 
[Kl]  terms  accounts  for  the  nonlinearities  arising  from 
large  displacements  and  interaction  between  those  displace¬ 
ments  in  the  reference  surface  and  the  small  strain  terms. 

Finally,  inclusion  of  the  [ K  2  )  terms  accounts  for 
nonlinearities  arising  from  interactions  taking  place  bet¬ 
ween  the  in-plane  reference  surface  deformations  and  rota¬ 

tions,  the  t h r ou gh - t h e - t h  i  c k n e s s  deformations  and  rotations, 
and  the  small  strain  deformations. 

In  practice,  inclusion  of  one  or  all  of  the  terms  (  K.0  ) 
through  ( K  2  J  is  a  simple  matter  and  is  left  to  the  discre- 

Some  examples  of  their  inclusion  are 


t 1  on  of  the  user. 


given  in  Chapter  IV. 


The  determination  of  { D ] ^  deserves  some  attention  as  the 


orientation  of  the  fibers  may  vary  from  layer  to  layer. 

Consider  a  single  orthotropic  layer  in  the  fiber  oriented 

reference  system  as  shown  in  Fig  (  5  ).  If  j  t*ien 


one  may  write: 
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1 1 


and  [ 0  ^  1  - [ C ) 


system 


t  o 


To  transfer  from  the  fiber  oriented 
the  global  system  requires  a 


tensor 


transformation  : 
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Following  this  transformation,  [D  )  is  seen  to  be 
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[  D  ^  ^  is  the  constitutive  matrix  in  the  local  fiber-oriented 

reference  system,  and  (D  |  is  the  matrix  in  the  global 

g  J 

system.  N'ote  that  ifii.  =0  or  90,  then  D16=D26“D36=D45=0. 


To  illustrate  the  formulation  of  <F>  we  return  to  the 
two  layer  example  and  note  in  detail  how  the  various  terms 
enter  in.  Considering  only  the  kth  node: 
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where  f  , f  , f  are  the  surface  tractions  in  each  of  the  u,v, 
u  v  w 


and  w  directions. 


If  one  has  on  1 y  f 


then  the  force  vector 


would  appear  as: 


F  ,  F^  , F  ^  , F  Q  are  then  the  generalized  forces  applied  to 
the  u.  f  f  degrees  of  freedom  at  the  k  t  h  node.  The 
limits  of  integration  p  are  either  or  depending  on  which 
face  the  load  is  applied. 

The  solution  of  the  strictly  linear  problem  was  a 
straight  forward  Crout  variation  on  the  Gauss  elimination 
method.  The  method  consists  of  a  factorization  of  the 
stiffness  matrix  into  the  product  of  a  lower  triangular 
matrix  and  an  upper  triangular  matrix.  The  use  of  this 
scheme  with  active  column  profile  (skyline)  storage  of  the 
stiffness  matrix  led  to  a  very  compact  program  and  allowed 
for  inclusion  of  a  resolve  capability  (i.e.,  new  load  cases) 
without  any  significant  additional  programming  effort.  (56) 


The  solution  of  the  nonlinear  problem  was  accomplished 

by  implementing  the  Newton-Raphson  or  modified  Niewton- 

Raphson  method  on  the  system  of  equations.  The  algorithm 
for  these  methods  is  essentially  the  same  with  the  exception 

that  the  stiffness  matrix  is  not  recalculated  at  each 

iteration  for  the  modified  method.  A  detailed  description 
of  these  methods  is  found  in  Zienkewicz  (41).  Appendix  III 
also  contains  further  discussion  on  the  nonlinear  solution 
technique  as  related  to  specific  test  problems. 

Chapter  III  is  devoted  to  an  explanation  in  further 

detail  as  to  how  the  computer  program  was  written  to  account 
for  the  integrations,  assemblies  and  solutions  of  the 

preceding  expressions. 


CHAPTER  III 


THE  FINITE  ELEMENT  ANALYSIS  PROGRAM 


To  validate  the  theory  presented,  a  finite  element 
analysis  program  was  written  to  solve  some  flat  plate 
problems  which  could  be  compared  with  previous  works. 

As  was  stated  in  the  theoretical  derivation  given  in 
Chapter  II,  every  effort  was  made  to  keep  the  computer 
program  element-independent.  That  is,  the  program  was 
written  to  use  any  element  in  the  family  of  2-D 
isoparametric  elements  ranging  from  the  simple  3-noded 
triangles  to  an  8-noded  Serendipity  quadrilateral  a  r.  d 
including  transition  elements  as  needed.  The  code  thus 
written  lends  itself  to  great  flexibility  in  element 
select  ion. 


Shape  function  routines  were  written  to  allow  the  user 
to  discretize  problems  quickly  and  reliably.  The  shape 
function  subprograms  evaluate  not  only  the  shape  function, 
but  also  its  derivatives  with  respect  to  the  global  coordi¬ 
nate  frame  (see  Fig.  a).  The  basis  for  the  subprograms  is 


the  4-noded  isoparametric  quadrilateral  where: 

■'■&)(  I  )/ 4- 


(102) 


and 


coordinates  of  the  nodes.  Making 


use  of  the  isoparametric  concept  yields  (adopting  the 


summation  convention): 


(103) 


(104) 


(  105  ) 


where  J  is  the  Jacobian  determinant  and  (  ), 

partial  derivative.  These  relations  define 
shape  function  routine  where  it  is  assumed 

coordinates  have  been  transformed  to  the 
coordinate  system. 


x  denotes  the 
steps  for  the 
that  the  nodal 
local  element 


The  extension  of  the  4-node  routine  to  higher  (or  lower) 


order  elements  is  accomplished  simply  by  adjusting  the  shape 


functions  to  account  for  the  addition  (or  deletion)  of 
nodes.  In  the  case  of  the  triangle  for  instance,  the  third 
and  fourth  nodal  shape  functions  are  added  together.  In  the 
case  of  the  8-noded  serendipity  quadrilateral  element,  each 
of  the  corner  nodal  shape  functions  is  adjusted  for  the 
presence  of  the  midside  nodes.  The  elements  which  one  may 
wish  to  use  as  a  transition  from  linear  to  quadratic  are 
obtained  simply  by  omitting  the  midside  node  for  the  linear 
edge. 


The  generation  of  the  stiffness  matrix  [K]  deserves 
special  attention,  since  many  zeros  exist  in  [M],  [L],  and 
[ D ] .  Concentrating  attention  on  a  single  node  i,  we  define 

[ B ]  = [ L |  [ M ]•  and  consider  the  matrix  triple  product 

T 

[K]^  =(B]. (D](B]j.  One  can  readily  see  from  the  dimensions 

of  the  matrices  that  formal  matrix  multiplication  would 
2 

require  2 4 M ~  + 1 9 2 M+ 3 3 0  multiplications  per  layer  where  M  is 

the  number  of  layers.  If,  on  the  other  hand,  one  were  to 

form  the  product  explicitly  accounting  for  the  zeros,  he 

would  find  that  the  number  of  multiplications  would  be 
2 

4M"+40M+7 I  per  layer.  For  instance,  given  M»3,  we  see  that 
one  would  perform  1122  multiplications  per  layer  formally  as 
compared  with  227  when  the  zeros  of  the  matrices  are  taken 
into  account.  Since  there  is  significant  difference  in  these 


numbers,  the  triple  product  was  formed  explicitly  and  then 


programmed . 


The 

choice  for 

the 

nodal 

degree 

of  freedom  vector  <a> 

a  s 

de  f i ned 

by  eqn  (SO) 

o  f 

C  h  a  p  t  e 

r  II  was  made  because  this  form 

becomes 

very  effic 

i  e  n  t 

when 

adding 

nodes  as  one's  choice 

o  f 

elements 

changes . 

For 

instance,  suppose  one  wishes  to  use 

a 

triangular  element. 

Then 

rilateral  element 
augmented  by  the 
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hand,  one  were  to 
,  the  <  a  >  and  {  Si  ] 
addition  of  a  and 


(106) 


choose  the  4-noded  quad- 
matrices  would  simply  be 
N  . 


Another  advantage  of  this  choice  for  (a>  has  been 
alluded  to  previously  in  the  way  that  [K]  is  formed.  Look¬ 
ing  at  the  triangular  element,  one  can  see  that  the  system 
of  equations  would  appear  as: 


(107) 


(108) 


The  program  was  written  to  take  advantage  of  this  submatrix 
structure.  Thus  the  whole  [K]  matrix  is  formed  by  succes¬ 
sive  calculation  of  ( K ] .  through  subscript  incrementation 
and  proper  stacking  of  each  submatrix.  Furthermore,  since 
[KJ  is  always  symmetric,  only  the  upper  triangular  portion 
need  be  calculated.  The  lower  half  was  formed  by  reflec¬ 
tion. 

In  order  to  facilitate  the  flexibility  of  varying  ele¬ 
ments,  it  was  necessary  to  incorporate  a  numerical  integra¬ 
tion  scheme  for  performing  the  necessary  integrations.  The 
Gauss-Legendre  quadrature  formulae  were  utilized  for  this 
purpose  since  "they  give  the  highest  accuracy  for  effort 
expended."  (48) 

This  choice  necessitated  that  the  limits  of  integration 
in  each  of  the  three  coordinate  directions  be  on  the  inter¬ 
val  from  -1  to  +1.  The  use  of  the  2-D  isoparametric  formu¬ 
lation  readily  accounted  for  the  in-plane  directions.  The 
t h r o ugh - t h e - t h i c kn e s s  coordinate  was  transformed  as  follows: 


(l  -  -1.) 
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where 

layer 

where 

Cion 


It  indicate 
,  and  /). 

the  layers 
becomes : 


s  the  layer,  t  is  the  thickness  of  the  kth 
is  the  transformed  variable.  In  the  case 

are  of  the  same  thickness,  the  transforma- 


z--  K  '  +  2t^L) 


and  thus 


ji  =  4-  ^ 
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(in) 


The  use  of  this  new  variable  yields  the  following  integ¬ 


ral  for  the  evaluation  of  [Kj 


i  J 
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Cast  in  terms  of  the  G a u s s - L e g e n d r e  quadrature  relations  one 
obtains : 


^  K-  r  5  t 


3  ) 


where  Q  ,Q  ,  (D  are  the 
r  s  t 

points 

and  fj  directions  and  p 
written  to  use  the  same 


weights  associated  with  the  sampling 
The  order  of  integration  is  n  in  the^ 
in  the  z  direction.  The  program  was 
integration  order  for  both  in-plane 


directions. 


,'rders  of  integration  for  those  directions  may 


CHAPTER  IV 


APPLICABILITY  AND  RANGE 

In  this  chapter,  the  applicability  and  range  of  the 
theory  are  presented.  The  first  application  examined  is 
that  of  a  plate  strip  (see  Fig.  7).  A  plate  strip  is 
infinitely  long  in  one  direction  (in  this  case,  the  y - 
direction).  The  simply  supported  strip  is  loaded  laterally 
(as  shown  in  Fig.  7)  so  as  to  cause  cylindrical  bending  to 
occur.  This  cylindrical  bending  implies  that  the  displace¬ 
ments  (w)  in  the  z -direction  are  symmetric  about  the  center 
line  of  the  strip  and  are  independent  of  y  and  displace¬ 
ments  (v)  in  the  y-direction  are  zero. 

The  first  case  studied  was  that  of  the  isotropic  plate 
strip  loaded  uniformly.  This  problem  posed  some  difficulty 
for  Witt's  (30)  parabolic  theory  using  the  constant  strain 
triangle  to  represent  the  reference  surface  displacements. 
The  decision  to  solve  this  problem  was  motivated  by  Witt’s 
difficulty.  With  the  present  cubic  theory  and  a  family  of 
finite  elements  to  choose  from,  a  study  could  be  carried 
out  to  investigate  how  Witt's  difficulties  could  be  over¬ 


come  . 


SECTION  A-A 


The  plan  of  action  was  to  see  how  the  cubic  theory 
compared  with  the  parabolic  theory.  This  was  accomplished 
by  solving  the  thin  isotropic  plate  strip  problem  using  the 
constant  strain  triangle.  Improvements  on  the  results  were 
expected  when  the  4-noded  and  8-noded  q u ad r i 1  a t e r a  1  s  were 
employed,  since  these  elements  would  permit  more  bending. 
The  logic  of  going  to  a  higher  order  interpolation  for  the 
in-plane  quantities  was  suggested  in  Witt's  (30)  work.  He 
noted  that  the  transverse  shear  strains  were  functions  of 
the  derivatives  of  the  in-plane  displacements  with  respect 
to  the  z  coordinate  and  the  derivatives  of  the  transverse 
displacement  with  respect  to  the  x  and  v  coordinates. 
Since  u  and  v  varied  quadraticallv  with  respect  to  z,  their 
derivatives  were  linear.  Since  w  varied  linearly  with 
respect  to  x  and  v,  that  derivative  was  always  a  constant. 
Therefore,  he  argued  that  no  matter  how  small  he  made  his 
element,  the  transverse  strains  would  never  go  to  zero,  and 
thus  there  would  always  be  this  shear  penalty  as  the  plate 
got  thinner. 

The  material  properties  for  the  plate  strip  were: 

E  ■  3.0  x  107  psi 
V  -  0.23 
h  *  0.2  inches 
a  ■  10  inches 


The  analytical  solution  for  this  problem 


was 


deter¬ 


mined  from  Fourier  analysis  using  Xavier's  method 


maximum  (w)  was  readily  found  to  be:  (55) 
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where  D  =  Eh  /  (  l  2  (  1  — ~  )  ] 


A  two-layered  finite  element  convergence  study  was 
performed  to  observe  the  way  the  elements  responded  to  the 


given  load .  Fig. 

8 

shows 

the 

various 

finite 

element 

arrangements  used. 

Note 

that 

the 

symmetry 

of  the 

problem 

was  exploited  to  limit 

the 

size 

o  f 

the  probl 

em.  Note,  also 

that  in  each  case  considered,  the  aspect  ratio  of  the 
elements  was  kept  equal  to  unity.  This  was  done  by  keeping 
the  dimension  (b)  in  the  y-direction  equal  to  the  length  of 
each  element  in  the  x-direction.  Such  a  technique  was 
possible  because  of  the  "strip"  nature  of  the  problem. 


5  8 


The  3-noded  (constant  strain),  -4-noded  (bi-linear),  8- 
noded  ( b  i - q u a d r a t  i  c  )  elements  were  studied  in  this  fashion. 
In  each  case,  two  sampling  points  were  used  to  perform  the 
integrations  in  the  z-direction  of  each  layer  (see  Fig.  6). 
Since  there  were  two  layers,  this  meant  integrations 
through  the  thickness  were  performed  with  -»  sampling 
points.  The  sampling  in  the  (^)  and  (^ )  directions  were 
varied  to  observe  the  reactions  of  the  elements  to  dif¬ 
ferent  integration  schemes.  For  the  constant  strain  and 
bilinear  elements  1  x  1  and  2x2  schemes  were  employed, 
while  for  the  quadratic  element,  2x2  and  3x3  schemes 
were  used.  Figs.  9-11  present  the  results  of  this  con¬ 
vergence  study. 

In  Fig.  9,  one  sees  the  results  for  the  constant 
strain  triangle.  Note  that  the  present  theory,  when  used 
in  conjunction  with  the  2x2  integration  scheme,  virtually 
duplicates  the  results  obtained  by  Witt.  The  only 


difference  is  that  in  Witts  work,  two  pseudolayers  were 
used  to  model  each  of  the  two  lamina,  where  the  present 
method  uses  no  pseudolayers.  I  n  other  words,  the  present 
theory  using  two  lamina  appears  equivalent  to  Witt's  method 
using  4  lamina.  Note  that  the  lxl  point  integration  scheme 
converges  more  quickly  early  on.  This  is  to  be  expected  as 
the  lxl  scheme  represents  a  "reduced"  (  -» 8  )  integration  for 
the  constant  strain  triangle.  The  2x2  scheme,  on  the  other 
hand  represents  an  "exact"  integration  method.  Both  inte¬ 
gration  schemes  suffer  from  the  same  drawbacks  as  did 
Witt's  3-noded  element.  That  is,  the  finite  element  solu¬ 
tion  converges  very  slowly  towards  the  analytical  results. 
The  best  solution  for  this  case  was  8  2  '1  of  the  analytical 
solution  obtained  with  80  elements  using  the  lxl  integra¬ 
tion  scheme.  Although  both  of  the  schemes  would  continue 
to  converge  to  the  exact  solution,  the  rate  of  convergence 
becomes  slower  and  slower  and,  therefore,  the  study  was 
terminated. 

Fig.  10  shows  the  comparisons  obtained  for  the  bi¬ 
linear  element.  Note,  again,  that  the  lxl  integration 
scheme  appears  to  converge  to  the  thin  plate  solution  very 
quickly  early  on  when  compared  to  the  2x2  scheme.  Again, 
as  with  the  triangle,  this  i s  to  be  expected  as  the  lxl 
scheme  is  a  reduced  integration  method.  (48)  There  was  a 
marked  improvement  in  the  convergence  of  this  element  over 


the  constant  strain  element  and  Witt's  triangular  element. 
The  best  solution  acquired  was  9  2%  of  the  analytical  results, 
obtained  using  90  elements  with  1  x  1  integration.  A  direct 
comparison  with  the  triangular  element  can  be  made  because 
the  same  number  of  nodal  points  (  and  thus  degrees  of  free¬ 
dom)  were  used  in  both  cases.  Again,  however,  the  conver¬ 
gence  was  very  slow,  and  a  prohibitive  number  of  elements 
would  be  required  to  obtain  100%  of  the  analytical  solution. 

The  results  obtained  with  the  2x2  integration  scheme  on 
the  s-noded  element  are  virtually  identical  to  those  obtained 
with  the  2x2  integration  on  the  triangular  element.  This 
leads  to  the  conclusion  that  any  improvement  hoped  for  by 
employing  the  quadrilateral  element  would  only  be  realized 
when  used  in  conjunction  with  the  lxl  integration  scheme. 

The  results  for  the  8-noded  biquadratic  element  are 
displayed  in  Fig.  11.  The  best  results  in  this  case  were 
obtained  by  using  2  x  2  integration.  Again,  about  92%  was 
the  best  this  element  could  do  with  the  isotropic  problem. 
The  3x3  integration  scheme  fell  far  short  of  approaching 
either  the  2x2  scheme  or  the  analytical  solution.  This 
method  appeared  to  level  out  at  about  74%  of  the  analytical 
solution  —  18%  below  the  2  x  2  method.  The  3  x  3  scheme 
represents  a  full  integration  method  which  has  been  shown  by 
other  authors  (for  instance  Cook  (36))  to  cause  this 


11  Convergence  of  8-node  Elements  —  Isotropic  Plate 


difficulty  as  the  plate  thickness  is  reduced. 


Since  the  2x2  integration  scheme  is  a  reduced  method 
of  integration  for  the  8-noded  quadratic  element,  the  re¬ 
sults  observed  in  Fig.  11  are  not  surprising.  When  the 
reduced  integration  scheme  of  the  8-noded  element  is  com¬ 
pared  with  the  reduced  scheme  for  the  -i-noded  element,  one 
sees  th3t  they  are  nearly  identical.  This  leads  to  the 
conclusion  that  improvements  in  performance  in  solving  the 
isotropic  problem  are  a  result  of  the  reduced  integration 
and  not  necessarily  a  result  of  the  higher  order  interpola¬ 
tion  function  for  the  reference  surface  displacements. 

It  should  be  stressed,  however,  that  in  each  case 
where  Witt  used  pseud Olivers  the  present  method  used  none. 
This  means  comparable  results  were  obtained  in  each  case 
using  full  integration  while  using  fewer  degrees  of  freedom 
than  Witt.  In  the  cases  where  reduced  integration  was 
used,  results  better  than  Witt's  were  obtained.  Also, 
since  only  one  fourth  as  many  integration  points  were  used 
for  the  reduced  scheme,  those  solutions  were  obtained  with 
about  one  fourth  the  computer  time  of  the  full  integration 
s  c  he  me  . 

In  every  case,  an  attempt  was  made  to  obtain  a  better 
solution  bv  increasing  the  number  of  lavers.  That  is, 


pseudolavers  were  tried.  The  use  of  these  pseudolayers  is 
comparable  to  increasing  the  number  of  elements  in  the  z 
direction.  In  every  case,  increasing  the  number  of  layers 
made  neglibile  improvement. 

A  review  of  the  literature  (51,52)  shows  that  this 
failure  of  the  "thick-plate"  element  to  converge  to  the 
"thin-plate"  solution  is  common.  These  "thick  plate"  theo¬ 
ries  are  essentially  of  the  Mind  1  i n  type  where  the  normal 
to  a  reference  surface  before  deformation  remains  straight 
but  not  necessarily  normal  after  deformation.  These 
assumptions  lead  to  stiffness  matrices  which  characterize 
bending  independent  of  the  transverse  shear  effects.  The 
difficulty  encountered  by  these  elements  as  the  plate  be¬ 
comes  thin  is  that  the  shear  stiffness  matrix  remains  too 
stiff.  Normally,  this  is  looked  upon  as  a  "penalty"  and  a 
number  of  techniques  are  formulated  which  attempt  to  reduce 
the  penalty.  The  methods  of  penalty  reduction  involve 
reduced  or  selective  integration  or  formulation  of  a  penal¬ 
ty  funct  .on  which  is  essentially  a  multiplier  on  the  shear 
stiffness  matrix  which  approaches  zero  as  the  plate 
thickness  is  reduced.  The  failure  of  the  thick  plate 
element  to  converge  to  the  thin  plate  solution  has  been 
attributed  to  the  "locking"  phenomena. 

The  present  formulation  is  not  of  the  Mindlin  type 


since  reference  surface  normals  are  allowed  Co  deform  after 


loading.  Furthermore,  the  derivation  of  the  stiffness 
matrix  is  carried  out  from  an  elasticity  point  of  view. 
The  stiffness  matrix  is  not  characterized  by  a  bending  part 
and  a  shear  part. 

The  thin-plate  solution  is  a  result  of  Kirchhorf's 
hypotheses  which  neglect  both  the  direct  strain  and  the 
shear  strains  >  tftj -2  '  In  addition,  the  normal  stress 
is  assumed  small  compared  to  J|(  and  so  that  it  is 
neglected  in  the  stress-strain  relations.  These 
assumptions  lead  to  a  theory  which  is  based  on  only  the 
transverse  displacement  w.  In  the  present  work,  the  direct 
transverse  strain  is  neglected,  but  the  shear  strains  are 
not.  The  problem  is  not  reduced  as  is  the  classical  thin 
plate  problem.  The  effects  of  the  transverse  shear  terms 
are  always  present. 

An  attempt  was  made  to  find  a  method  for  reducing 
the  effects  brought  about  by  the  shear  strains.  The  method 
arrived  at  was  to  lower  the  magnitudes  of  the  shear  moduli  1 


Gn* 

and 

G  2  3  ’ 

This  reduction  resu 

lted  in  a  so 

ftening  of 

the 

element  and  a 

convergence  to  the 

thin-plate 

solution. 

Fig- 

1  2 

shows  the 

effects  of  varying 

the  reduction  factor. 

The 

vert 

ical  axis 

shows  the  choices 

for  factors 

to  use  to 

multiply  G  ^  ^  and  G1^.  The  horizontal  axis  shows  the  resul- 


tant  maximum  w  displacement  for  the  isotropic  strip  prob¬ 
lem.  Note  that  the  thin  isotropic  plate  solution  (w=l)  is 
obtained  when  a  multiplication  factor  of  1/100  was  used. 

Looking  back  at  the  "locking"  problem  and  methods 
which  have  been  used  to  eliminate  the  locking,  one  sees 
that  in  each  case  there  was  an  attempt  to  soften  the  ele¬ 
ment  in  some  way  so  as  to  eliminate  the  effects  of 
transverse  shear.  In  effect,  what  was  being  done  was  to 
artificially  reduce  the  effects  of  G  1  3  and  G  2  3  on  the 
solution  of  the  problem.  The  present  work  takes  a  less 
artificial  and  more  direct  approach  to  the  problem.  Since 
the  element  was  designed  for  a  thick  structure,  the  three 
dimensional  effects  had  to  be  considered.  Direct 
transverse  stress  and  transverse  shear  strains  were  always 
considered.  The  only  way  to  reduce  their  effects  was 
through  the  constitutive  relation  (D). 


Since 

this 

case  was  that  of 

a  thin  isotropic 

plate 

a  /  h  =  5  0  )  , 

i  t 

lent  itself  well 

to  observing  the 

way  a 

normal  line  segment  would  deform  under  the  load.  A  thin 
plate  line  element  normal  to  the  reference  surface  should 
come  ciose  to  duplicating  the  assumptions  of  Kirchof f  '  s 
hypothesis.  That  is,  the  normal  should  remain  straight  and 
perpendicular  to  the  reference  surface.  fig*  13  is  a 
plot  of  u  vs.  z,  which  shows  that  Kirchhoff's  hypothesis  is 
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W  =  1  REPRESENTS  THE 
THIN  PLATE  SOLUTION 


NOTE  THAT  THE  THIN  PLATE 
SOLUTION  IS  OBTAINED  WHEN 
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13  Check  of  Klrchhoff's  Hypothesis  —  Isotropic  Plate 


being  met  by  the  cubic  assumption  for  u  in  the  z-direction 


The  rotations  (  ^ • >  are  such  that  the  normals  do  indeed 

remain  straight  and  normal  after  loading. 


The  second  case  studied  was  identical  to  the  first. 


except  that 

orthotropi 

c 

materials  were 

layers.  The 

plies  were 

oriented  at  0 

uniformly  at 

163.84  ps 

i  . 

The 

materia 

were: 

Ei 

2.5  x 

1  0  '  p  s  i 

E2 

= 

1.0  x 

6 

10  p  s  i 

<*12> 

3 

0.25 

G12 

3 

5.0  x 

6 

10  p  s  i 

G  1  3 

= 

5.0  x 

,  .  6 

1  U  p  S  1 

G2  3 

3 

2.0  x 

106  p  s i 

h  *  0.2  inches 

a  ■  10.0  inches 


The  analytical  solution  for  this  problem  was  deter¬ 


mined  as  before.  The  only  difference  being  that  D“Q  h  / 1 2 


where  Q 


1  1 


1  ,  2 
:  / [ e  -  V  ei. 

1  1  1  12  2  1 


Figs.  14-16  show  the 
the  elements  of  interest, 
to  the  exact  solution  was 
each  case,  the  reduced 


convergence  study  carried  out  for 

In  nearly  each  case,  convergence 
obtained.  Xote,  also,  that  in 

integration  scheme  converged  more 
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rapidly  than  the  full  integration  scheme.  In  the  cases  of 
the  3-noded  and  4-noded  elements  using  full  integration, 
the  convergence  to  the  thin  plate  solution  was  obtained 
only  when  the  dimensions  of  the  elements  were  reduced  to 
about  one  third  of  the  thickness  of  the  element.  This  is 
totally  unacceptable  from  a  practical  standpoint  and  led  to 
the  conclusion  that  for  this  thin  plate  problem  reduced 
integration  of  the  3-noded  or  4-noded  element  was  impera¬ 
tive. 


The 

only  case 

which 

did  not  converge  was 

that  of  the 

q  uad  r  a  t i c 

element 

using 

3x3  integration. 

A  s 

with  the 

isotropic 

problem, 

the  e 

lement  again  appeared 

t  o 

converge 

to  a  so 

lution  wh 

ich  was  nearly  20%  below 

the 

solution 

obtained  by  the  2x2  scheme  and  the  exact  solution.  The 
3x3  integration  scheme  appears  to  make  the  8-noded  element 
too  stiff.  This  problem  seems  to  be  common  with  the  8-noded 
element  and  is  referred  to  specifically  in  reference  (56). 

The  choice  of  element  is  always  left  to  the  discretion 
of  the  user.  It  appears,  however,  that  the  8-noded  element 
with  2X2  integration  is  the  most  desirable  based  on  this 
convergence  study.  Even  a  single  element  produced  results 
which  were  99%  of  the  exact  solution. 


As  with  the  isotropic  case 


u  vs  z  was  plotted  at  x»0 


in  order  to  observe  the  deformation  of  the  normal.  Fig-  17 
shows  that  the  normal  remained  undeformed  and  normal  to  the 
reference  surface.  This  was  in  accord  with  the  expecta¬ 
tions  of  the  author  since  the  plate  was  thin  ( S  =  5  0  )  and  the 
plies  were  both  oriented  at  0  degrees. 

The  third  case  studied  was  again  cylindrical  bending 
of  a  strip  plate.  This  time,  however,  different  ply  lay¬ 
ups  were  studied  in  order  to  observe  the  element's  perform¬ 
ance  for  symmetric  as  well  as  asymmetric  problems.  Compar¬ 
isons  were  made  with  analytical  results  obtained  by  P3gano 
(  7  )  .  More  specifically,  three  ply  lay-ups  were  examined  — 
(0,0),  (0,90),  and  (0,90,0).  The  material  properties  used 
in  the  previous  case  were  also  employed  here.  For  this 
case,  however,  the  loading  was  changed  from  uniform  to 
sinusoidal.  That  is, 

P  (  x  )  5  Pq  s  i  n  (if  x  /  a  )  (lti) 
Pg=162.7io  psi 

The  reason  for  choosing  sinusoidal  load  was  to  check 
the  solution  against  Pagano's  (7)  solution.  Additionally, 
this  provided  a  check  on  the  ability  of  the  program  to 
handle  a  more  complex  loading  situation. 

The  finite  element  arrangement  shown  in  Fig.  8(i)  was 


used  for  this  case. 


That  is, 


five 


8-noded  elements  were 


used,  and  2X2  integration  was  employed.  The  boundary  con- 


di t ions 

were  the  same 

as  in  previous 

cases. 

The 

element  1  s 

response  to  the 

problem  of 

plates  of 

various 

thicknesses 

was  of  interest. 

In  order 

that  the 

reader 

may  better  vi 

sualize  the  plat 

es  and  gain 

a  perspec- 

t  i  v  e  of 

"thickness"  , 

refer  to  Fig. 

18.  This 

represents 

values 

for 

S=a/h  from 

50  to  4.  S  =■  1  0  is 

considered 

(for 

0  r  t  h  0  t 

r  0  p  i  c 

plates)  to 

be  a  thick  plate, 

because  the 

ef  - 

f  e  c  t  s 

0  f 

shear  stress 

result  in  displacements  which 

are 

much  greater  than  predicted  by  the  classical  laminated 
plate  theory. 

Figs.  19--21  are  plots  of  w"  vs.  S  for  the  various  ply 

—  3  14 

orientations,  where  w  =  (  l  0  0  E  _  h  /  (  p  .  a  )]w  .  These  figures 

compare  the  finite  element  results  with  those  obtained  from 
Pagano's  elasticity  solution,  a  finite  element  solution  by 
Owen  (32)  and  the  classical  laminated  plate  theory  (CLPT). 
In  each  instance,  the  finite  element  method  duplicated  the 
elasticity  results  showing  that  the  element  has  application 

over  a  wide  range  of  thicknesses  for  orthotropic  plate 
problems.  The  figures  also  serve  to  amplify  the  statement 

made  in  Chapter  I  that  the  transverse  shear  deformations  must 
be  considered  when  analyzing  thick  composite  plates. 
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To  gain  a  perspective  of  why  the  transverse  shear 
cie  formation  becomes  more  and  more  important  as  plate  thick¬ 
ness  is  increased,  it  is  instructive  to  observe  Figs.  22 — 
24  which  present  plots  of  TT  vs.  "z  at  x  =  0  ,  where  “'E^u/hp^ 
and  ? = z / h . 

Fig.  22  displays  a  very  thick  (S*4)  asymmetrically 
oriented  (0,90)  plate.  Note  that  the  F.E.  method  follows 
very  closely  the  results  obtained  by  Pagano  through  the 
elasticity  approach.  Figs.  23  and  24  show  u  vs.  z  for  the 
symmetric  ply  lay-up  (0,90,0).  In  each  case,  the  finite 
element  method  follows  the  elasticity  solution  closely. 
The  influence  on  thickness  is  readily  observed  by  compari¬ 
son  of  these  two  curves.  In  the  first  (  S  ■  1  0  )  ,  the  normal 
deformation  effect  due  to  shear  is  seen  to  be  quite  close 
to  the  CLPT  solution.  As  the  thickness  increased  to  S  =  >«  , 
the  normal  deformation  became  much  m  o  r  ■»  apparent  and  one 
readily  sees  how  this  deformation  becomes  more  important  as 
the  thickness  is  increased. 

Figs.  23  and  24  also  show  the  cubic  nature  of  the 
approximation  for  u(z).  Note  particularly  the  inflection 
point  in  the  middle  ply  at  z*0.S.  The  "smooth"  nature  of 
the  approximation  is  also  readily  observed  in  these  fig¬ 
ures.  The  interlaminar  boundary  discontinuity  of  Pagano’s 
solution  is  not  duplicated  because  of  the  cubic  approxima- 


8w 


t  i  o  n  . 


Th  i  s 


smoothing  has  an  impact  on  the  evaluation 


u,  and  thus  on  (£  )  near  the  boundary  (as  shown  in  Figs 
z  ***■ 


2  6  —  2  8)  . 


Figs.  25  —  2  8  show  plots  of  (£j^)  vs.  T  at  x»0  where 

/p  .  In  the  case  of  a  thick.  (S“4)  two-ply  plate 
U 

with  fibers  oriented  at  (0,0),  one  sees  (Fig.  25)  that  the 
F.E.  values  follow  quite  closely  those  pre  ticted  by  the 

elasticity  solution.  When  the  fiber  orientation  is  changed 
to  (0,90),  there  is  a  small  degradation  in  the  results  for 
the  shear  near  the  interlaminar  boundary  (see  Fig.  26).  For 
the  3-ply  (0,90,0)  cases.  Figs.  27  —  28  show  that  the  F.E. 
method  follows  the  exact  solution  well,  except  at  the 
boundaries.  In  Pagano's  solution,  (  ,  u  ,  a  n  d  w 

were  kept  continuous  across  the  interlaminar  boundaries. 
In  this  F.E.  solution,  the  continuity  is  imposed  on  u  and 
u  with  w  constant  through  the  thickness.  These  differen¬ 
ces  in  boundary  conditions  at  the  interfaces  account  for 
the  differences  between  the  two  solutions.  In  other  words, 
by  imposing  an  assumption  for  the  continuity  of  the  first 
derivative  of  u  with  respect  to  z  at  the  interlaminar 
boundries,  the  theory  is  seen  to  violate  equilibrium  at 
those  discrete  points  through  the  thickness  where  abrupt 
changes  in  material  properties  occur.  In  terms  of  the 
stress  resultants  which  are  obtained  by  integration  of  the 
stresses  through  the  thickness,  equilibrium  is  satisfied. 


CLPT  — 


PAGAN'O 


It  should  be  noted  that  these  stress  results  were 
obtained  at  the  Gauss  integration  points  of  the  1st  element 
and  extrapolated  linearly  to  the  x  »  0  boundary.  This  extra¬ 
polation  also  may  contribute  to  some  of  the  differences 
observed  •  .  the  interlaminar  boundaries. 

Figs.  29  —  33  show  plots  of  ( )  vs.  "z  at  x  =  a/  2  where 
(  Cfx  )  =  (‘fjt  )  /  p  ^  .  In  each  case,  the  elasticity  solution  is 
followed  nearly  identically. 

The  fourth  case  to  be  presented  is  that  of  the  square 
orthotropic  plate  shown  in  Fig.  3a.  The  plate  was  modelled 
using  25  8-noded  elements  for  the  quarter  plate.  Reduced 
(2x2)  integration  was  employed.  The  material  properties 
were  again  the  same  as  in  the  previous  orthotropic  cases. 
Loading  in  this  problem  was  a  transverse  sinusoidal  load  of 
the  form: 

p  =  p^sin(Jfx/a)sin(ry/a)  (116) 

The  boundary  conditions  for  this  plate  were  simple 
supports  all  around.  Along  the  edges,  the  transverse  dis¬ 
placement,  w,  was  fixed,  and  the  normals  to  the  reference 
surface  at  the  edges  were  not  allowed  to  rotate  perpendic¬ 
ular  to  the  boundary. The  u  and  v  displacements  were  free. 
The  plate  was  a  symmetric  (0,90,90,0)  lay  up. 
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Fig.  35  shows  a  comparison  of  w  vs  S  for  this  plate. 
The  finite  element  solution  compares  nearly  identically 
with  the  elasticity  solution  and  gives  results  which  are 
better  than  another  F.E.  solution  by  Pryor  and  Barker  (25). 
Thus,  the  element  appears  to  respond  well  to  changes  in 
geometric  restraints  and  configurations. 

To  test  the  nonlinear  theory,  the  symmetric  (0,0)  lay¬ 
up  and  the  asymetric  lay-up  (0,90)  cases  of  the  plate  strip 
problem  were  examined.  Since  the  present  theory  allows  the 
inclusion  of  two  sets  of  nonlinear  terms  (Kl  and  K2),  a 
check  was  made  to  determine  the  influence  each  of  these 
terms  had  on  the  solution.  As  with  the  linear  case,  5  8  - 
noded  elements  (with  2x2  integration)  were  used  to  model 
the  plate  strip. 

Fig.  36  shows  plots  of  v7  vs  Q  for  both  the  (0,0)  and 
the  (0,90)  cases.  The  length  to  thickness  ratio  for  these 
cases  was  chosen  to  be  S“10,  which  means  that  this  is  a 
thick  plate  strip.  As  expected,  the  (0,90)  strip  is  more 
flexible  than  the  (0,0)  strip.  The  solid  line  depicts  the 
solution  which  was  obtained  by  including  only  the  [  K  0  ]  or 
linear  terms.  The  lines  represented  by  the  squares  show 
the  results  of  including  the  [K1J  terms  which  contain 
nonlinearity  with  respect  to  the  reference  surface.  The 
inclusion  of  the  f K2 ]  terms,  which  contain  n o n  1  i  n e a r  i  t  l  e s 


with  respect  to  the  thickness,  is  depicted  by  the  triangle 
symbols  on  the  lines.  Vote  that  when  the  [ K 1 ]  terms  are 
included,  there  is  a  stiffening  of  the  plate  strip.  That 
is,  there  is  less  displacement  for  a  given  load.  Mote  also 
that  when  the  [ K  2  ]  terms  are  included,  the  solution  is 
stiffer  than  the  linear  solution  but  more  flexible  than  the 
solution  obtained  using  the  ( K l )  terms. 


In  Figs.  19  and  20  it  was  shown  that  when  [ K  0  )  terms 
were  used  in  the  analysis,  the  solutions  obtained  were 
identical  with  the  classical  laminated  plate  theory.  Fur¬ 
ther,  it  was  shown  that  as  the  plate  became  thicker,  the 
presence  of  the  shear  terms  in  the  thickness  direction 
produced  transverse  displacements  (w)  which  were  greater 
than  the  classical  laminated  plate  theory.  Since  the  [KO] 
terms  contain  linear  expressions  for  u,v,  and  w  (the  refer¬ 
ence  surface  displacements)  and  linear  expressions  for  the 


rotations  < 


$jx  > 


ich  account  for  the  effects  of  the  shear 


terms,  [KO]  was  able  to  effectively  represent  the  elastici¬ 
ty  solution.  This  leads  to  the  conclusion  that  when  shear 
terms  are  taken  into  account,  there  is  an  accompanying 
increase  in  the  transverse  displacements.  The  inclusion  of 
{  K  1  )  in  the  analysis  takes  into  account  nonlinear  displace¬ 
ments  in  the  reference  surface  but  the  terms  for  the  ef¬ 


fects  of  transverse  shear  are  linear  as  they  were  in  (KO). 


Because 


of  the  presence  of  these  nonlinear  terms  for  the 


reference  surface  displacements,  the  solution  becomes  s  t  l  t  - 


fer  in  comparison  with  the  solution  obtained  using  only 
(  K  0  ]  .  When  [  K  2  ]  is  included  and  the  fully  nonlinear  formu¬ 
lation  is  used,  the  "softening"  phenomena  occurs.  The 
softening  is  a  direct  result  of  the  important  nonlinear 
rotation  terms  which  are  contained  in  (K2). 

Figs.  17  and  38  show  the  errects  of  including  (FI]  and 


(  K 2  ]  on 

the  f  and 

*  X  x  z 

s  t  r  e  s 

ses. 

That  is,  the  stresses 

increase 

when  the  [  K  2  ] 

terms 

are 

included. 

The 

sixth  case  pre 

s  e  n  t  e  d 

1  s 

again  a  square  ortho  tropic 

plate 

(see  Fig.  3s). 

The  material  p 

ropert les 

w  e  r 

e  re 

t  a  i  n  e  d 

from 

previous  cases. 

The  load  for 

this  case 

wa 

s  d 

i  s  t  ri¬ 

b  u  t  e  d 

uniformly  in  the 

transverse  di 

rection  as 

P  =  P0 

The  boundaries  of  the 

plate  were  all 

clamped  ( 

r  1 

mp^d 

me  a  n  s 

that 

all  degrees  of  f 

reedom  at  the 

boundaries 

w  e 

r  e 

f  i  x  e  d  ) 

and  the  piv  lay-up  was  (0,0).  The  plate  was  modelled  by 
using  16  s-noded  elements  (lxl  integration)  for  one  quad¬ 
rant  instead  of  the  8-noded  element  with  2x2  integration 
simply  to  conserve  computer  time.  This  problem  was  chosen 
as  a  check  of  how  the  non-linear  terms  of  the  theory 
would  respond  to  a  change  in  geometry.  The  check  was  made 
by  proportionally  increasing  the  load  p  ^  and  observing  the 
transverse  deflection  (w)  at  the  center  of  the  plate  where 
it  was  a  maximum.  Fig.  39  shows  a  graph  of  the  results 
obtained.  This  plot  gives  a  comparison  between  the  present 


work., 


the  work  of  Witt  ('10),  and  .in  approximate  analytical 
solution  obtained  by  Chi  a  (10).  Vote  in  particular  that 
the  inclusion  of  the  [  K  1  ]  terms  gives  results  which  are 
some wh. it  better  than  those  obtained  by  Witt  who  included 
[Kl]  and  nearly  identical  to  Chia’s  (who  used  a 
perturbation  method  on  the  Von  Karman  plate  equations). The 
inclusion  of  the  higher  order  terms  of  [ K  2  ]  yields  an 
expected  softening  of  the  element  as  the  loads  were  in¬ 
creased.  The  same  arguments  as  were  put  forth  for  the 
plate  strip  problems  apply  here. 
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This  work  represents  the  first  time  that  the  inclusion 
of  these  (  K  2  ]  terms  has  been  performed  and  therefore  no 
direct  comparison  to  previous  work  can  be  made. 


For  a  detailed  account  of  computer  time  and  technique 
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used  in  solving  these  nonlinear  problems  refer  to 


Appendix 


CHAPTER  V 


COSCLL'SIOS’S 

A  number  or  conclusions  have  been  made  in  the  Ust 

chapter  in  conjunction  with  the  discussions  of  the  results. 

One  of  the  significant  contributions  of  this  work,  lies 
in  the  unique  use  of  cubic  spline  theory  which  has  been 

well  known  for  many  years  in  the  field  of  interpolation  and 

in  the  solution  of  differential  e  <>u  ations.  The  most  recent 
work  with  cubic  splines  in  these  fields  has  always  concen¬ 
trated  on  the  displacements  <u>  being  known  and  solving  for 

the  rotations  {$)  in  terms  of  the  known  displacements. 

I* 

Here  (u>  is  referred  to,  because  of  the  context  of  the 

present  use  of  spline  theory.  In  reality,  no  one  to  the 

best  of  the  author's  knowledge  has  heretofore  considered 
the  cubic  spline  to  model  through-the-thickness  displace¬ 
ments.  In  all  previous  uses  of  the  cubic  spline  in 
conjunction  with  plate  solutions,  the  spline  was  used  to 

model  in-plane  displacements.  In  this  work,  the  cubic 

spline  equations  are  solved  in  the  reverse  sense.  That  is, 
the  <u>  is  treated  as  the  unknown  and  are  solved  for  in 

terms  of  the  (  £  > .  This  process  resulted  in  the  develop¬ 

ment  of  a  cubic  shape  function  which  did  not  introduce  any 
new  degrees  of  freedom  over  the  previously  derived  para- 


ho  1  i  c  representation.  The  use  of  the  spline  in  this 
"reverse"  role  is  a  promising  new  technique  in  numerical 
approximation  theorv  where  the  rotations  ( become  the 
degrees  of  freedom  and  <u>  is  then  represented  in  terms  of 
them  through  the  derived  shape  function. 


A  consequence  of  the  use  of  an  odd-degree'd  spline 


(cubic  in  this  case)  is  that  the  equations  necessary  to 
represent  the  interpolation  functions  are  very  compact  and 
easily  represented.  A  look  at  Palazotto's  work  (57), 
wherein  he  used  (essentially)  a  parabolic  spline,  shows  the 
complexity  of  the  expressions  which  results.  In  terms  or 
computer  implementation,  this  is  important  not  only  in 
terms  of  computational  effort,  but  also  in  terms  of 


programmer  effort.  The  cubic  spline,  although  a  higher 
order  approximation,  leads  to  a  simplification  of  the  equa¬ 
tions  and  lends  itself  to  less  programmer  effort.  The 
reduction  in  programmer  involvement  results  in  fewer  er- 


The  results  clearly  show  that  the  use  of  the  cubic 
spline  theory  to  represent  t h r o u g h - t h e - t h i c k n e s s  displace¬ 
ments  in  thick  composite  plates  is  a  viable  and  effective 
application.  Comparisons  with  previous  works  show  that 
there  is  the  ability  to  obtain  equally  accurate  results 
with  fewer  degrees  of  freedom.  Consider  for  instance  work 


I 


or  Epstein  (42)  (where  a  number  of  pseudolayers  was  neces¬ 
sary  to  obtain  the  results  presented).  This  same  technique 
was  used  bv  Witt  in  his  work.  With  the  cubic  spline,  each 
layer  can  be  represented  as  a  layer,  and  not  a  number  of 
pseudolayers. 


A 

comp 
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use 

of  3-D  e 

each  1 

aver 

again  clearly  shows 

the  redo 

of  d  e 

g  r  e  e  s 

of  freedom  and 

t  h  u 

s  comput 

time. 

I  n 

both  the  present 

work  and  Pa 

there 

are 

2(M+l)+3  degrees 

o  f 

f  reedosi 

tlie  number  of  layers).  In  Epsteins  (->2)  formulation  there 
are  5  (  M  +  1  )  degrees  of  freedom  per  node.  With  a  variety  of 
ID  brick  elements  to  choose  from  there  is  a  range  from 
1  ( M  + 1  )  to  5 ( M+ 1  )  degrees  of  freedom  per  node  (or  even  more) 
depending  on  the  choice  of  brick  element  used.  When  pseu¬ 
dolayers  are  used,  there  are  essentially  2M  layers  used  to 
represent  M  layers.  In  numerical  terms,  it  one  were  to 


consider 

a  6-lavered 

problem 

,  one  would 

find 

that 

Epstein'  s 

element 

with  pseud 

olayers 

won  Id  r  e  q  u  i 

re  39 

DOF 

per  node. 

Witt's  e 

lement  would 

require 

29  DOF  per 

node, 

the 

tril inear 

brick,  element  with  no  pseudolayers  would  require  21  DOE  per 
node,  and  the  present  element  would  require  17  DOE  per 
node. 


The  difficulties  in  converging  to  the  thin 


isotropic 


plate  solution  were  not  overcome  by  using  a  higher  order 
element  for  the  in-plane  interpolations.  Reduced  integra¬ 
tions  of  the  bi-linear  quadrilateral  and  the  bi-quadratic 
quadrilateral  improved  the  results  from  82 'i  to  92'.  of  the 
thin  plate  solution.  The  thin  isotropic  plate  solution  was 
obtained  when  the  effects  of  the  transverse  shear  terms 
were  reduced  by  reducing  the  transverse  shear  moduli  i  by  a 
factor  of  100. 

The  inclusion  of  the  [ K  2  )  terms  in  the  nonlinear 
formulation  is  an  advance  in  the  state  of  the  art  in 
through-the  thickness  representations.  It  is  seen  to  become 
more  important  as  load  levels  are  increased  (assuming  yield 
stress  is  not  exceeded).  Since  [ K  2  ]  contains  nonlinear 
transverse  rotational  terms,  it  also  would  become  more 
important  whenever  transverse  shear  terms  affect  the  prob¬ 
lem.  The  thicker  the  structure,  the  more  important  the 
inclusion  of  these  terms  becomes. 

The  inclusion  of  the  [ K  2  )  terms  leads  to  somewhat 
larger  transverse  stresses  than  those  obtained  with  (  K  1  ]  . 
This  is  due  to  the  increase  in  the  transverse  displacements 
which  accompanies  the  inclusion  of  the  [ K  2 ]  terms. 


In  summary,  the 

cubic 

assumpt ion 

for  the  variations 

o  f 

the 

/ 

displacements  u 

and  v 

through  the 

thickness  proved 

t  o 

be  an  effective  approach.  In  terms  of  degrees  of  freedom, 
there  was  actually  a  decrease  over  the  parabolic  assumption 
simply  because  accurate  results  could  be  obtained  without 
resorting  to  pseudolayers.  Secondly,  difficulties  of  lock¬ 
ing  when  applying  this  thick  plate  element  to  a  thin  plate 
were  overcome  by  reducing  the  values  for  G  ^  and  G1 ^ . 
Thirdly,  the  continuously  differentiable  nature  of  the 
cubic  assumption  for  u  and  v  through  the  thickness  resulted 
in  discontinuities  in  the  shear  stresses  at  the  interlami¬ 
nar  boundaries.  This  is  seen  as  a  limiting  factor  on  the 
applicability  of  the  element.  It  should  be  noted,  however, 
that  if  the  stresses  were  integrated  through  the  thickness 
and  resultant  forces  were  considered  then  equillibrium  is 
satisfied  within  the  element.  Fourthly,  the  inclusion  of 
the  nonlinear  terms  of  [  K  2  ]  which  represent  nonlinear  ef¬ 
fects  arising  through  the  thickness  leads  to  greater  dis¬ 
placements  at  a  given  load  level  than  the  inclusion  of  just 
the  reference  surface  nonlinear  terms  of  [  K 1  ]  .  In  all,  it 
is  felt  that  the  elements  developed  (namely  the  -t-noded 
with  lxl  integration  and  the  8-noded  with  2x2  integration) 
are  viable  and  accurate  when  compared  to  both  the  analyti¬ 
cal  solutions  and  other  finite  element  solutions. 
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APPENDIX  I 

Based  on  the  development  of  the  nonlinear  strains  in 
Chapter  II,  the  following  expressions  show  the  results  of 
making  substitutions  for  the  hVs  and  the  "el  ^'s  and  then 
simplifying  to  obtain  the  expressions  for  the  Lagrangian 
strains  appearing  in  the  boxes. 
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The  expressions  Tor  [LO] ,  [Ll] ,  and  (L 2]  are  p 
e  for  completeness: 
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APPENDIX  III 


This  appendix  is  devoted  to  a  more  detailed  d  e  s r  i  p  - 


t  i  o  n 

o  f 

the  element  as 

well  as 
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more 
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rough 
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the 

computer  times 

and  tec 
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lques 

use 

d  to 

solve 

the 

nonlinear  problems. 

As  discussed  in  the  body  of  the  text,  the  finite 
element  can  be  considered  to  be  a  n  u  a  s  l  -  1 D  element.  That 
is,  it  consists  of  i  2-1)  isoparametric  element  containing 
degrees  of  freedom  of  u  ,  v ,  and  w  at  each  of  its  nodes. 
This  2-D  element  is  used  to  disc  retire  the  reference  sur- 
rioe.  The  effects  of  the  deformations  through  the  thick¬ 
ness  are  accounted  for  by  including  rotational  degrees  of 
freedom  at  the  reference  surface,  the  outer  surface  of  the 
structure,  and  at  each  of  the  interlaminar  boundaries. 
This  leads  to  an  element  which  has  2  M  +  5  degrees  of  freedom 
per  node  as  can  be  seen  in  Fig.  -*2. 
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is  cussed  in  Chapter  IV  were  solved 
Ea  ch  of  the  nonlinear  problems 
using  the  Newton-Raphson  method, 
by  an  incremental  load  method 
phson  iteration  between  each  load 
turn.  This  incremental-iterative 
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The  solution 
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were  generated  by  fir 

s  t 

solving 

the 

linear  problem,  i 
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using  { K  0  J  terms  onl 

y  • 

The  load 

was 

then  incremented 

by 

the  amount  which  is  i 

n  d  i  c 

a  t  e d  by 

the 

graphical  symbols 

o  f 

Figs.  36  and  39.  For 

Pagano's  problem 

using  a  (0,0)  ply 

lay 

-up,  the  incremental 

load 

factor 

was 

25.  For  the  same 

pro 

blem  but  with  a  (0,90) 

ply 

lay-up  , 

the 

incremental  load  factor  was  chosen  to  be  10.  The 
incremental  load  factor  used  in  Chia's  problem  was  chosen 
to  be  50.  These  choices  were  made  based  on  observation  of 
the  linear  solution  at  the  maximum  load  level  desired  and 
selected  so  as  to  give  favorable  graphical  points  for  the 
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To  illustrate  the  rate  of  convergence  of  the  solutions 
for  each  of  the  load  levels.  Table  III.l  is  presented. 
This  table  shows  how  the  solution  progressed  for  Pagano's 
(0,0)  ply  lay-up  problem.  Note  that  the  inclusion  of  the 
full  nonlinear  stiffness  matrix  terms  [  K  2  ]  did  not  result 
in  a  drastic  increase  in  the  number  of  iterations  required 
to  reach  equi 1  librium  at  each  load  level.  There  was  only  a 
nominal  increase.  These  solutions  were  obtained  in  nearly 
the  same  computer  time  as  well.  The  [ K 1 ]  solution  was 
obtained  in  1085  CP  seconds  while  the  [K2]  solution  was 
obtained  in  1107  CP  seconds. 
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TABLE  I  I  I  .  1 


Comparison  of  Iterations  Required  for  Convergence 


Load  Level 

Using  [ K  1  ] 

Using  [ K  2  ] 
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1  5  0 

2  1 

2  2 

1  7  5 

2  4 

2  5 

2  00 

2  6 

2  7 

The  convergence  criteria  which  was  used  to  determine 
equilibrium  was  suggested  by  Zienkiewicz  (-*8).  Iteration 
was  continued  until  the  norm  of  the  force  imbalance  vector 
was  less  than  some  tolerance  times  the  maximum  observed 
norm  tor  previous  iterations.  For  the  nonlinear  problems 
shown  in  this  work,  the  tolerance  was  set  at  .0001. 

Table  III. 2  shows  the  progression  of  the  solution  for 
Chia's  problem.  Again,  there  was  no  large  increase  in 
computer  resources  required  to  include  the  [K2]  terms.  The 
[K1J  solution  (which  was  solved  to  a  load  level  of  200)  was 
obtained  using  62  CP  seconds  while  the  [ K  2  )  solution  (which 
was  solved  to  a  load  level  of  300)  was  obtained  using  318 


CP  seconds. 
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abstract 

A  non-linear  thick  composite  shell  theory  is  presented 
in  which  the  through-the-thickness  displacements  are  model¬ 
led  ns  .  ng  a  variation  of  a  cubic  spline.  The  theory  is 
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restricted  to  small  s  t  r  uni.  The  imposition  of  the  cubic 
distribution  thro  ugh- the -thick ness  insures  that  the 
comparability  of  the  displacements  and  their  first  and 
second  derivatives  and  thus  the  shear  strains  are  maintained 
from  lamina  to  lamina.  The  cubic  distribution  is  seen  as  a 
higher  order  approximation  than  has  been  previously 
employed,  but  because  of  the  nature  of  the  spline,  the 
theory  is  less  cumbersome  and  more  easily  implemented  than 
the  parabolic  theory.  In  addition,  there  is  no  introduction 
of  additional  degrees  of  freedom  with  the  cubic  theory.  A 
family  of  2-D  i s o p a r a m e t r  i  c  elements  is  employed  in 
conjunction  with  the  theory  to  solve  a  class  of  3-D  thick 
plate  problems.  Results  are  presented  showing  comparisons 
which  are  in  good  agreement  with  previous  work. 
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